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Impact of neuronal heterogeneity on correlated colored 
noise-induced synchronization 
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Synchronization plays an important role in neural signal processing and transmission. 
Many hypotheses have been proposed to explain the origin of neural synchronization. 
In recent years, correlated noise-induced synchronization has received support from 
many theoretical and experimental studies. However, many of these prior studies have 
assumed that neurons have identical biophysical properties and that their inputs are well 
modeled by white noise. In this context, we use colored noise to induce synchronization 
between oscillators with heterogeneity in both phase-response curves and frequencies. 
In the low noise limit, we derive novel analytical theory showing that the time constant 
of colored noise influences correlated noise-induced synchronization and that oscillator 
heterogeneity can limit synchronization. Surprisingly, however, heterogeneous oscillators 
may synchronize better than homogeneous oscillators given low input correlations. We 
also find resonance of oscillator synchronization to colored noise inputs when firing 
frequencies diverge. Collectively, these results prove robust for both relatively high noise 
regimes and when applied to biophysically realistic spiking neuron models, and further 
match experimental recordings from acute brain slices. 
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1. INTRODUCTION 

Synchronization of neural oscillators is thought to play a criti- 
cal role in sensory, motor, and cognitive processes (Sanes and 
Donoghue, 1993; Fries et al., 2001; Wang, 2010). In many net- 
works, synchronization is achieved via direct coupling such as 
through gap junctions and chemical synapses. However, there 
are several systems (notably, the mammalian olfactory bulb) 
where the mode of coupling is less clear and neural synchrony is 
hypothesized to arise from partially correlated presynaptic inputs 
(Galan et al., 2006; Marella and Ermentrout, 2010). Indeed, in 
non-oscillatory networks of neurons, such correlated input is 
largely responsible for the output correlations of the neurons 
(de la Rocha et al., 2007). Thus, a natural question is: how do the 
properties of neurons and networks alter output correlations for 
a given degree of input correlation? At small input correlations, 
output and input correlations can be regarded as linearly pro- 
portional; this ratio is called the susceptibility (Shea-Brown et al., 
2008). For example, (de la Rocha et al., 2007) showed that the 
susceptibility depends on the background firing rate of the neu- 
ron. For some model systems, this susceptibility can be computed 
using linear response theory (which assumes small perturbations 
around the stationary state). 

When neurons fire regularly, they can be regarded as noisy 
nonlinear oscillators and, as such, there are many mathemat- 
ical techniques available for their analysis. In particular, the 
phase-response curve (PRC) provides a compact and useful char- 
acterization of the responses of a nonlinear oscillator to external 



perturbations. The PRC describes the shift in timing of, say, an 
action potential as a function of the timing of the input rel- 
ative to the last action potential. In several studies, we have 
described the theoretical relationship between the shape of the 
PRC and the ability of identical neurons to transfer partially syn- 
chronized activity (Marella and Ermentrout, 2008; Abouzeid and 
Ermentrout, 2009). In these studies, the only source of hetero- 
geneity considered between neural oscillators was their unshared 
(uncorrelated) inputs, which consisted of white noise. Recently, 
we extended these methods to cases in the low noise limit in 
which the oscillators were not identical and showed how het- 
erogeneity in intrinsic properties could significantly degrade the 
output correlation in pairs receiving common inputs (Burton 
et al, 2012). 

In this study, we extend this theory to include colored noise 
inputs and, further, report some surprising effects of heterogene- 
ity. First, we derive a set of equations for the distribution of 
phase differences for pairs of heterogeneous oscillators driven by a 
partially correlated Ornstein-Uhlenbeck (OU) process (low-pass 
filtered noise). We next show that the theory developed for phase 
reduced models works well with a conductance-based biophys- 
ical model. We then show that, quite surprisingly, at low input 
correlations, heterogeneity can sometimes produce higher out- 
put correlations than the homogeneous case. That is, consider 
two distinct oscillators, A and B, such that the AA pair has a 
small susceptibility and the BB pair a larger susceptibility. Then, at 
low correlations, the susceptibility of the AB pair can sometimes 
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exceed that of the AA pair. We confirm this somewhat counter- 
intuitive prediction with recordings from regularly firing mitral 
cells of the main olfactory bulb. In addition to heterogeneity in 
response properties, neurons can fire at different frequencies, and 
such frequency differences can significantly impact correlated- 
noise induced synchronization (Markowitz et al, 2008; Burton 
et al., 2012). Here, we find that for some frequency differences 
between oscillators, there is an optimal time scale of correlated 
noise that will maximally synchronize the oscillators. We do not 
see this effect when the oscillators have the same frequency. 

2. MATERIALS AND METHODS 

2.1. PHASE REDUCTION MODEL 

In Appendix, we provide a brief overview of how to reduce a gen- 
eral weakly perturbed limit cycle to a single differential equation 
for the phase of the cycle. If we assume that the original limit cycle 
represents repetitive firing of a single compartment neuron model 
that is driven by a noisy current, I(t), then we obtain: 

^ = 1 + 6A(6)/(t)/C m (1) 
at 

where C m is the membrane capacitance, 6 is the phase (or, typ- 
ically, the time since the last spike), and A (6) is the PRC of the 
neuron. The PRC describes the phase-dependent shift in the spike 
times of an oscillator receiving small perturbations. It is readily 
measured in neurons and other biological oscillators (Torben- 
Nielsen et al., 2010) and provides a compact representation of 
the effects of stimuli on the timing of action potentials. A (9) 
has dimensions of milliseconds per millivolt; that is, the shift in 
timing of the next action potential per millivolt perturbation of 
the potential. Mathematically, for a given model, A (6) is found 
by solving a certain differential equation (see Appendix). It is a 
periodic function of phase and, with no loss in generality, we can 
normalize the period to be 2jt for simplicity. 

2.2. STATIONARY DENSITY 

Given the reduced model (Equation 1), we can now turn to the 
main question at hand, which is: how do oscillating heteroge- 
neous neurons transfer correlations? We will consider two types 
of heterogeneity: differences in the PRC shapes and differences 
in natural frequencies. We drive the oscillators with correlated 
filtered noise. After reduction to phase variables, we obtain: 



e; = 1 +€Ai(6i)x (2) 

0' 2 = l+tA 2 (e 2 )y + € 2 W (3) 

x' = -x/x + % X /Jt (4) 

/ = -y/x + $ r /Vt (5) 



01 and 62 are the phases of two oscillators, and Ai(9) and A 2 (9) 
are PRCs of the two oscillators. Without loss of generality, we set 
the natural frequency of one oscillator to 1. The parameter to then 
determines the magnitude of the difference in natural frequencies 
between the two oscillators. 6 <gC 1, thus the noise is weak and the 
frequency difference is small. The processes x and y are generated 



by an OU process with the same time constant x. % x and %y are two 
correlated white noise processes, i.e., (^ x (t)% x (t')} = 8(f — f')> 
<^(f)^(f')> = 8(f - O. (%x(t)%y(t')) = c8(f - t'), where c is the 
degree of correlation. 

We remark that the allowable frequency difference is 0(c- 2 ), 
which seems considerably smaller than the magnitude of the 
noise, which is t. However, as the noise has zero mean, what mat- 
ters is the variance of the noise, which has magnitude c- 2 . Thus, 
the scales of both the frequency difference and the synchroniz- 
ing inputs (correlations in the noise) are similar. If the frequency 
differences are larger, then no synchronization is possible. 

Our goal is to compute the stationary distribution of the 
phase difference between two neurons since this will enable us to 
compute various measures of correlation and synchrony. Thus, 
some variable substitution will be helpful: 9 = 9i, cp = 9 2 — 9i. 
Therefore, cp is the phase difference between the two oscillators. 
With this change of variables, the equations are: 

0' = 1 + €Ai(9)x (6) 
(t>' = f[A 2 (0 + (t))y- Ai(0)x] + € 2 w (7) 

and x, y are as above. Let p(x, y, 9, cp, t) represent the probability 
density function at time t : 

Pr(X(t) &{x,x + dx), Y(t) &(y,y + dy), ©(f) € (9, 9 + dQ), 
$(f) e (c(>, <|> -f- deb)) = p(x, y, 9, §)dxdyd$d§ (8) 

We denote the stationary density (long-time behavior asf^- 00) 
as p ss (x, G, cp). 

Our goal is to compute the probability density of the phase 
difference between the two oscillators, R(<\>), which is: 

/ / p S5 (x, y, 9, c|>) dxdydQ (9) 

-00 J— 00 J 0 

If the oscillators are perfectly synchronized, then R($) will be a 
delta function centered at cp = 0. If the oscillators are completely 
independent, then -R(cp) = l/(2n). In Appendix, we show that 
R($) satisfies a simple first order boundary value problem (BVP). 
We present the exact equation for this in Results. 

2.3. ORDER PARAMETER 

Once we get the distribution of phase differences, R(<\>), we need 
a number to estimate the synchronization, which means the 
sharpness of this distribution. In this study, we use an order 
parameter (OP) to do this. We define: 

op = Vc 2 + S 2 (10) 

p2tz 

C= R(<$>) cos(cp)d(t) 
Jo 

S= / J?(cp) sin(<p)c?cp 
Jo 

9 = atan2(C, S) 
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OP is a representation of sharpness and 6 is the estima- 
tion of the peak position. For certain types of heterogeneity, 
R(<&) is peaked at <\> = 0; in this case, we can show that the 
cross correlation of the spike times is (R(0) — l/(27t))/(27t) 
(Burton et al., 2012). However, OP provides a better global 
measure of the synchrony and is not dependent on the peak 
being centered at 0; we will therefore use OP in our current 
results. 



for the model is the voltage component of the solution to 
the adjoint equation (Equation 32). The software package XPP 
(Ermentrout, 2002) includes an algorithm for computing the 
adjoint solution for an exponentially stable limit cycle, so we 
simply compute various limit cycles (say with very different 
parameters but similar periods) and then compute Ami (6) for 
those specific parameters. We save the result as a lookup table and 
then numerically solve the phase equation. 



2.4. M0RRIS-LECAR MODEL 

The Morris-Lecar (ML) model (Rinzel and Ermentrout, 1989) 
is a simplified two-dimensional system membrane model that 
we use to compare the phase models with a full biophysical model: 



= h - gdV x - V L ) - gxWiiVi - V K ) 

at 

- gCa fWoo(Vi)(Vi - V Ca )+OX 

dwi _ Woo(Vi) - Wi 
~dt~^ Tw(Vl) 
dV 2 

C-— =h- gL(V 2 - V L ) - g K w 2 (V 2 - V K ) 
at 

- gc a mcc(V 2 )(V 2 - V Ca ) + ax 
dw 2 _ Woo(V 2 ) - w 2 

~dt~^ x w (V 2 ) 
x' = - x/x + % x /*Jx 



(11) 
(12) 

(13) 

(14) 

(15) 
(16) 



with <^(t), %i{t')) = Ht-f), (%i(t), i= 2 (f0> = Ht-f), and 
(%i(t)i ^2(f')> = c Ht — t'), c G [0, 1]. The auxiliary functions are: 



moo (V) = 0.5 • (1 + tanh((V - V a )/V b )) 
W(X (V) =0.5 • (1 + tanh((V - V c )/V d )) 

r w (V) 1 



cosh((V- V c )/(2V d )) 



(17) 
(18) 

(19) 



84 mV, V L = 

mS 



The parameters used in this paper are: Vk 
-60 mV, V Ca = 120 mV, g K = 8^, g L = 2*%, g Ca = ±— 2 , 
C = 20^, V a = -1.2 mV, V b = 18 mV, V c = 2 mV, and V d = 
30 mV. Ii, I 2 and 4>i , § 2 vary for each figure. 

To get the phase from the noisy voltage signal generated by 
the ML model, we first apply the Hilbert transform (HT) to V(t) 
which allows us to get a phase. However, the phase is not uni- 
form as it is not a temporal phase. We then map the HT phase to 
a temporal phase on the noise-free limit cycle which gives a uni- 
form phase-distribution as required by the theory. This allows us 
to estimate R($) for the biophysical model, where 4> here is the 
phase difference between two ML model neurons that are driven 
with partially correlated noise. 

In some of the figures, we simulate the phase-reduced dynam- 
ics for the ML model. To do this, we must compute the 
infinitesimal PRC, Ami (9)- As described in Appendix, the PRC 



2.5. NUMERICS 

To get solutions to the stochastic phase and membrane equations, 
we use the Euler-Murayama method. We solve the BVP for the 
stationary phase difference density using a custom BVP solver 
written in MATLAB. All codes are available by request. 

3. RESULTS 

3.1. APPROXIMATION OF THE PHASE DIFFERENCE DENSITY 

Oscillators driven with a correlated fluctuating signal will exhibit 
a degree of synchronization that depends on the size of the sig- 
nal, the strength of correlation, and the similarity of the two 
oscillators. Thus, for example, identical oscillators driven by 
small enough identical white noise will synchronize perfectly 
(Pikovsky et al, 1997; Teramae and Tanaka, 2004). The rate 
at which these identical oscillators synchronize depends on the 
properties of the noise - in particular, its autocorrelation (Nakao 
et al., 2007; Goldobin et al, 2010). In general, and especially 
in biological systems, there will be a great deal of heterogene- 
ity in any pair of oscillators. For example, for neurons, there 
is always some source of independent noise so that the input 
correlation is always less than 1. The neurons may also be fir- 
ing at slightly different frequencies. Finally, even if the neurons 
are adjusted to fire at the same frequency, their distribution of 
ion channels can be very different and, thus, their response to 
correlated signals can be quite different (Burton et al., 2012). 
If the fluctuating inputs are sufficiently small, then any sta- 
ble limit cycle oscillator can be reduced to a so-called phase 
model where the dynamics are characterized by a single vari- 
able, the phase, such that the firing is considered to occur at 
a phase of 0 and the time between spikes is mapped onto an 
angle between zero and 2tt. Here, we consider driven pairs of 
heterogeneous oscillators that receive partially correlated filtered 
noise. As our main examples come from neuroscience, we assume 
that the external inputs are implemented as currents, in which 
case the phase model for the pair of neural oscillators has the 
form: 

e; = l + eAiCeoxCf) 

6' 2 = l + € 2 w -MA 2 (e 2 )y(t) 

where x(t) and y(t) are OU processes with the same time 
constant, x, and with correlation c; Ai 2(6) are the PRCs for 
the two oscillators; € is a small positive parameter (charac- 
terizing the magnitude of the fluctuations); and go accounts 
for the frequency difference in the unperturbed oscillators (see 
Materials and Methods, Equations 2-5). We are primarily inter- 
ested in the distribution of the phase difference, $ := Q 2 — 
01. In the Appendix (Equation 62), we show that R(<$>), the 
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probability density function for the phase difference, satisfies a 
simple BVP: 



d<$> 



{[c • g($) - Ci]Rm + (4nw - C 2 )R(<$>) = K 



R($) = R($ + 2%) 
g(4>) = g(* + 2it) 
R($)d$ = 1 



K = 2u> 



2tt 



The 2ti— periodic function g(<\>) and the constants, C\ t 2, depend 
in a complicated way on the forms of the PRCs and the time con- 
stant of the noise, x (see Appendix). However, all quantities can 
be found by integrating elementary functions. If the oscillators 
have the same PRC, then C2 = 0 and g($) is even symmetric. 
If the oscillators have the same frequency, then co = 0. When 
both C2 and w vanish, we can immediately solve the BVP, yield- 
ing R($) = N/(C\ — cg((\>)), where N is a normalization constant 
so that the integral is 1. This is the result found in Marella and 
Ermentrout (2008) for white noise, but is clearly also true for 
colored noise. When the oscillators are identical and there is no 
difference in frequencies, the phase difference density is symmet- 
ric and always peaks at 0. However, when 00 — Ci is nonzero, 
the peak of the phase difference density will generally be off- 
set. We note that € does not appear in the expression for R($), 
which says that the phase difference density is, to a first approx- 
imation, independent of the amplitude of the noise. Figure 1 
shows typical results comparing the perturbation calculation and 
the simulation of the Langevin equation. In Figure 1A, at fairly 
high noise 6=1, there is some distortion at the peak of the 
distribution, but as predicted from the theory, the distribution 
magnitude is largely independent of the magnitude of the fluctua- 
tions. Figure IB shows a similar simulation, but the correlation of 
noise is lower (c = 0.5 vs. c = 0.8), the noise is faster (x = 0.25 vs. 




FIGURE 1 | Novel analytical theory of correlated colored noise-induced 
synchronization of heterogeneous oscillators matches Monte Carlo 
simulations for low to moderate levels of noise. Stationary phase 
difference density is shown as computed from the solution of the BVP and 
through Monte Carlo simulation from t = 1000 to f = 201000 in steps of 
0.05. Monte Carlo data binned into 100 bins between —it and u. There is a 
frequency difference of e 2 /2 where c is the magnitude of the noise. Here 
Ay(9 ; ) = sin(a,0 — sin(9i + at) + £>, sin(28 ( ), where / = 1 , 2 for two 



oscillators. (A) x = 1 , a q = 0.1, a 2 = 0.6, b\ = 0.32, 62 = 
(B) x = 0.25, at = 32 = 0.5, b\ —bi = 0.3, and c= 0.5. 



0.3, and c = 0.8. 



x = 1), and the PRCs are identical. In this case, even the higher 
noise simulations match the theory. We once again emphasize 
that the perturbation expansion requires a small value of e, but 
clearly, the simulations show that 6 can be nearly 1 and still yield 
good agreement. 

We note that the density of the phase differences can be related 
to more conventional measures of correlation. In Burton et al. 
(2012), we showed that the spike time cross-correlation (CC) 
between a pair of weakly noisy oscillators is: 



CC(t) 



1 

2jt 



R(-x) 



1 

2jt 



(20) 



For example, if the oscillators are asynchronous, then they have 
a uniform phase difference density and the cross-correlation 
will be 0. This calculation confirms ones intuition that differ- 
ent neurons that receive correlated noise will have spike time 
cross-correlations that peak off-center. 

Figure 2 shows that we can apply the theory through two 
levels of simplification. The ML system is a simple, biophysi- 
cally realistic model for a spiking neuron (Rinzel and Ermentrout, 
1989). With different choices of parameters, the onset to oscilla- 
tory behavior can be either through a Hopf bifurcation (HB) or a 
saddle-node on an invariant circle (SNIC) bifurcation. The PRCs 
that result from these distinct bifurcations are often quite differ- 
ent (Brown et al, 2004; Izhikevich, 2007) and thus have quite 
different synchronization properties. In Figure 2, we tune the ML 
model so that each cell has the same frequency but the parameters 
are quite different and so the PRCs are different (see Figure 2B). 
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FIGURE 2 | Analytical theory accurately predicts synchronization of 
biophysically realistic spiking neuron models. (A) Invariant phase 
difference density computed from the reconstructed phase of two ML 
model neurons receiving partially correlated colored noise (period is 
91.25 ms, x = 5 ms, c = 0.8). Three cases are illustrated with either 
identical (homogeneous) or mixed (heterogeneous) PRCs. The "Hopf" case 
corresponds to a set of parameters where the oscillatory activity arises via 
a HB and the "SNIC" case through a SNIC bifurcation (Rinzel and 
Ermentrout, 1989); see Appendix for parameters. (B) PRCs for the two 
cases. (C) Same as (A), but using simulations of the phase reduced 
equations. (D) Solutions to the BVP using the PRCs from (B). 
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In Figure 2C, we show the results of a Monte Carlo simulation in 
which the biophysical model is driven by correlated noise. Phase 
is reconstructed from the voltage traces using a Hilbert transform 
and from these, we obtain phase difference histograms. In this 
figure, the correlation c is 0.8, x = 5 ms, and the natural period 
of the oscillation is 91.25 ms. For the same degree of correla- 
tion, two HB oscillators are much better at synchronizing than 
are two SNIC oscillators. This result is consistent with the the- 
ory developed in Marella and Ermentrout (2008) for white noise 
and also for spike time correlations over fast timescales (i.e., spike 
synchrony) (Barreiro et al., 2010). At this high correlation, the 
heterogeneous HB-SNIC pair shows greatly reduced synchrony 
from either of the homogeneous cases and a shift in the peak 
even though there is no frequency difference. Figure 2B shows the 
two PRCs that were determined using the adjoint method. We 
then used these PRCs to compute the invariant densities for the 
corresponding phase reduced models. The invariant density is 
a function that describes the distribution of phase differences 
of the two neurons over some time interval consisting of many 
cycles. Thus, the peak of the invariant density indicates the most 
likely phase difference, and a large peak at zero phase differ- 
ence would indicate that the two neurons are well synchronized. 
Comparison between Figures 2A,C shows excellent agreement. 
Finally, we substituted the numerically computed PRCs into our 
BVP and computed the invariant density. The result of this cal- 
culation is shown in Figure 2D. There are small differences in the 
amplitude, but the shapes and the shift of the densities in the het- 
erogeneous case are almost identical. Thus, through two levels of 
reduction (first, from the full model to the phase model, and sec- 
ond, from the Langevin phase model to the approximate invariant 
density), we see that our analytical method works very well at esti- 
mating the invariant density of phase differences between neural 
oscillators. 

3.2. PRC HETEROGENEITY 

Our approximation of the invariant density, while requiring that 
we solve a BVP, allows us to explore the effects of heterogeneity 
much faster than simulating the appropriate Monte Carlo system. 
Thus, we will use this method to explore the effects of PRC het- 
erogeneity, frequency differences, and the color of the noise on the 
ability of oscillators to synchronize. One simple global measure 
of synchrony/correlation for systems whose natural dynamics 
are periodic is the circular variance, <3 c ir<k = 1 — OP, where we 
define an order parameter (OP) (see Materials and Methods, 
Equation 10): 



OP = 



£ 



R($) cos fydfy) + 



r 



R(§) sine)) d<$> 



For a flat distribution, OP = 0 (a c ii-cle = 1) and for a delta func- 
tion distribution, OP = 1 (a c ircle = 0). The OP is a commonly 
used measure for the degree of synchronization between two 
oscillators (Kuramoto, 2003). 

In general, one expects that the synchrony between two oscil- 
lators forced with correlated noise would be greatest if the oscil- 
lators are homogeneous. Certainly, if the inputs are identical 



(i.e., no independent or unshared noise), then identical oscillators 
will synchronize perfectly, while heterogeneous oscillators will not 
synchronize perfectly. That is, the phase density will not be a 
delta function. [See Burton et al. (2012) for a proof]. However, 
surprisingly, at low input correlations, it is possible for a hetero- 
geneous pair of oscillators to produce greater synchrony than one 
(but not both) of the respective homogeneous pairs of oscillators. 
Figure 3 illustrates the behavior of two separate homogeneous 
pairs of oscillators (blue and green lines, respectively) as the input 
correlation varies from 0 to LA third, heterogeneous pair com- 
prised of an oscillator from each homogeneous pair is shown in 
red. Figure 3A shows the two different PRCs; pairs of oscillators 
with the green PRC ("PRC 1-PRC 1") produce weaker synchrony 
than pairs of oscillators with the blue PRC ("PRC 2-PRC 2"). 
This is demonstrated in Figure 3B, where the correlation is set 
to 0.8. Note that the phase difference density for PRC 2-PRC 
2 pair is more peaked than that for PRC 1-PRC 1 pair, while 
both densities are more peaked than the heterogeneous "PRC 
1-PRC 2" pair. As noted above, the peak of the heterogeneous 
pair is not at the origin but rather, is shifted to the left. In order 
to get a global measure of synchrony, we plot OP as a function 
of the input correlation (Figure 3C). As c -¥ 1, both homoge- 
neous pairs approach OP = 1 (i.e., perfect synchrony) while the 
heterogeneous pair never exceeds OP = 0.4. However, at low cor- 
relations, the heterogeneous pair can actually synchronize better 
than the PRC 1-PRC 1 pair (compare red to green lines in inset). 
That is, in the presence of low correlations, a "good synchro- 
nizer" paired with a "bad synchronizer" performs better than the 
homogeneous pair of bad synchronizers. This effect is not just due 




FIGURE 3 | Oscillator heterogeneity can enhance correlated 
noise-induced synchronization at low input correlations. (A) Two PRCs 
with the form A/(8/) = sin(ay) - sin(a ; + 6,) + b/sin(20j), si =0.1, 
bi = 0.32, s 2 = 0.6, and b 6 = 0.3. (B) fl(ctj) with different combinations of 
PRCs. Blue: PRC 2-PRC 2; red: PRC 1-PRC 2; green: PRC 1-PRC 1. Solid 
lines are theoretical predictions from the solution to the BVP and open 
symbols are Monte Carlo simulation results (same notation applies in 
following figures). Parameters used: t = 1 and c= 0.8. (C) Synchronization 
as input correlation varies from 0 to 1; inset shows magnification when 
c < 0.5. (D) Same as (C), but using the ML model. Parameters used: 
/, = 110, §i = 0.04616, l 2 = 120, and §2 = 0-04. 
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to our approximate expansion as the Monte Carlo simulations 
show the same phenomenon. Figure 3D further hints that we can 
also see the effect in the full ML model, although the results are 
not as clear. 

3.2.1. Experimental evidence 

Could this subtle difference in the ability of neural oscillators to 
transfer correlation be seen in experiments? To answer this, we 
re-examined data from a previous study (Burton et al., 2012). 
Mitral cells from the mouse main olfactory bulb were injected 
with constant current overlaid with frozen noise to evoke noisy 
periodic firing. PRCs were then experimentally estimated using 
our previously described method using the spike-triggered aver- 
age (Ermentrout et al, 2007). [Complete methods are provided 



in Burton et al. (2012)]. In this dataset, we found several exam- 
ples where injecting partially correlated noise produced greater 
synchrony between two different mitral cells firing at the same 
rate than for one of the mitral cells across different trials (experi- 
mentally simulating a homogeneous pair of mitral cells). Figure 4 
illustrates an example. In Figure 4A, we show the voltage traces 
(top) of two mitral cells receiving correlated inputs, and the spike 
times (middle) and phase (bottom) as determined by a simple lin- 
ear interpolation between spikes. Figure 4B shows the PRCs from 
each of these two cells along with their fit to the exponential-sine 
PRC model (see Appendix, Equation 64). In Figure 4C, we show 
the phase difference density as constructed from the linear phase 
interpolation of the two cells. In this example, the currents deliv- 
ered through the electrodes are perfectly correlated. However, 
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FIGURE 4 | Physiological neuronal heterogeneity can enhance correlated 
noise-induced synchronization at low input correlations. (A) Example 
linear interpolation of phase between recorded spike times of two mitral cells 
injected with perfectly correlated colored noise. Top: experimentally recorded 
membrane potentials. Middle: raster plot of spike times. Bottom: phase. (B) 
Experimentally estimated PRCs for the two cells shown in (A). Dashed lines 
are fits of the exponential-sine PRC model to the estimated PRCs. (C) Phase 
difference densities of the two cells during injection of perfectly correlated 
currents. Densities were calculated from pairs of 5 sec recordings. Blue and 
green curves show densities for homogeneous pairs of cell 1 and cell 2, 
respectively. The red curve shows the density for the heterogeneous pair of 
cell 1 with cell 2. (Di) Experimental and (Dii) theoretical OP vs. input 
correlation for homogeneous and heterogeneous pairs of the two mitral cells. 
Theoretical curves calculated by solving the BVP with the model PRC fits and 
x = 5. [Note that the same results were obtained in separate calculations for 
x = 3, the time scale of the noise used in Burton et al. (2012)]. (Ei-Eii) Mean 
OP (±SEM) vs. input correlation across 85 pairs of mitral cells (formed from 
27 separate mitral cell recordings described in Burton et al. (2012)). For each 



pair, the cell with the greatest area under its homogeneous OP vs. correlation 
curve was classified as the "good synchronizer" of the pair. (Fi) Theoretical 
OP vs. input correlation (with x = 3) for each of the 85 heterogeneous pairs 
from (E) (light red lines), plotted against the theoretical OP vs. input 
correlation of a homogeneous pair formed from the average mitral cell PRC. 
Note that many (but not all) of the heterogeneous pairs exceed the 
homogeneous pair in the low correlation range shown. On average (dark red 
line), physiological heterogeneity enhances synchrony for input correlations 
up to ~0.27. (Fii) Magnification of the homogeneous and average 
heterogeneous lines from Fi for low input correlations. (Gi) Percent and (Gii) 
absolute change in theoretical OP for heterogeneous vs. homogeneous bad 
pairs of mitral cells. Black lines plot OP changes for pairs in which 
heterogeneity increased synchrony at low input correlations; magenta line 
plots mean OP enhancement (±SEM) for these pairs. Grey lines plot OP 
changes for pairs in which heterogeneity did not increase synchrony. Note 
that heterogeneity mediates the greatest percent increase in OP at low 
(<0.1) input correlations, similar to experimental results shown in (E). Inset 
shows magnification when |AOP| < 0.03. 
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unlike the simulations, the neurons themselves are intrinsically 
noisy, so there is a substantial component of "private" noise. 
Nevertheless, one can see that cell 1 (blue) synchronizes better 
across trials than does cell 2 (green) across trials. Figures 4Di,Dii 
show the OP as reconstructed from the experimental data and as 
obtained by using the computed PRCs, respectively. This shows 
that at low correlations, the heterogeneous pair ("1-2") can syn- 
chronize better than the "2-2" homogeneous pair (but not the 
"1-1" homogeneous pair). The inset in 4Dii magnifies the low c 
region. 

Are the results presented in Figures 4A-D for a single pair of 
mitral cells consistent across a larger population of mitral cells? 
To answer this, we examined recordings from 27 regularly fir- 
ing mitral cells, from which we were able to form 85 different 
pairs of mitral cells with highly similar (<5Hz difference) fir- 
ing rates. For each pair of mitral cells, we computed the OP 
across varying input correlations for both homogeneous com- 
binations and the heterogeneous combination. We automatically 
classified the mitral cell with the greatest homogeneous OP across 
all levels of input correlation as the "good synchronizer" of the 
mitral cell pair. Figure 4E shows the mean OP vs. correlation 
across the 85 good, bad, and heterogeneous mitral cell pairs. 
Note that, even with this relatively insensitive classification of 
good vs. bad synchronizers, there is a region at low input corre- 
lations where, on average, heterogeneous pairs synchronize better 
than the bad homogeneous pairs. This phenomenon is seen more 
clearly when we use the experimentally estimated PRCs and the 
BVP to compute the OP vs. input correlation. Figure 4Fi plots 
OP vs. c for all heterogeneous pairs (light red lines), the mean of 
the heterogeneous pairs (dark red line), and the OP for a single 
homogeneous pair whose PRC is the mean of all the PRCs (black 
line). For many cases (but not all), heterogeneity increases the OP 
above that achieved by a uniform population of neural oscilla- 
tors with the mean PRC. Figure 4Fii magnifies the mean OP vs. c 
curves at low correlation; the red curve is clearly higher than the 
black curve. 

We then quantified the degree to which physiological levels of 
heterogeneity [as experimentally measured in mitral cells (Burton 
et al., 2012)] can enhance synchrony between neural oscillators. 
Using the BVP and our experimentally estimated mitral cell PRCs, 
we calculated the percent and absolute change in OP for all 85 het- 
erogenous vs. homogeneous bad mitral cell pairs. That is, for the 
example pair in Figure 4Dii, we subtracted the green from the red 
line to calculate the absolute change in OP, and divided this dif- 
ference by the green line to calculate the percent change in OP. 
Figures 4Gi,Gii plot the results of this analysis for all 85 pairs. 
In 26 of these pairs (plotted in black), heterogeneity enhanced 
synchrony at low input correlations, with a mean increase in 
OP (plotted in magenta; ±SEM) of up to 36%. Thus, in rela- 
tive terms, physiological levels of heterogeneity can significantly 
enhance correlated noise-induced synchrony at low input correla- 
tions. While this relative enhancement in synchrony corresponds 
to an admittedly low absolute increase in OP of up to 0.01 on 
average (Figure 4Gii), we nevertheless expect this phenomenon 
to significantly contribute to patterns of oscillatory synchrony 
in the olfactory bulb and potentially other brain regions (see 
Discussion). 



3.2.2. Good vs. bad synchronizers 

When is a neuron a good vs. bad synchronizer? Here, the BVP is 
much simpler since we just have to compare homogeneous pairs. 
In this case, the probability density function can be written as: 



*(<!>) 



N 



-g(0) 



(21) 



where N is a normalization and g((\>) is defined above by setting 
n = m. For low values of c, we get: 
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Since the two neurons are identical, the peak of R(<$>) occurs at 
4> = 0 and, so, we can estimate the zero lag cross-correlation as 
[R(0) — 1/(2jt)]/(2ji). Using the approximations above, we see 
that: 

' 2 * g($) 
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That is, the cross-correlation is linearly proportional to the input 
correlation (for small c) and this factor [called the susceptibil- 
ity (de la Rocha et al., 2007)], is a simple function of g(<\>). We 
can maximize S if we can make the integral as small as possi- 
ble. Note that g($) is periodic, and the integral over a period 
is proportional to the constant Fourier coefficient. Recall that 
g(<\>) is a low-pass filtered version of h(§), which is the auto- 
correlation function of the PRC. Thus, h(0) is positive and so 
is g(0). The integral of g($) is proportional to the integral of 
fc(4>), which is just 2i\.a^ where an is the DC component of 
the PRC. Hence, we can minimize the integral and maximize 
the correlation transfer (susceptibility) if we mimimize the DC 
component of the PRC. This fact generalizes the conclusions in 
Marella and Ermentrout (2008) and Abouzeid and Ermentrout 
(2009) that state that more sinusoidal PRCs are the best syn- 
chronizers. For example, with a PRC of the form (sin(a) — 
sin(x + a)) exp(C(% — lit)), we obtain the best synchrony when 
a = — arctan C. 

Can we determine when a pair of oscillators will have the prop- 
erty that a good-bad heterogeneous pair is better than a bad-bad 
homogeneous pair? Since the effect is only seen at low corre- 
lations, this suggests a perturbation expansion for small c. We 
write -R(4>) = ^ c n R n ($) and find that Rq is constant and so to 
order 1, R(§) = Rq + cRi(<\>). From this, we can compute OP, 
OP = c L cos 4>.Ri (())) d<\>. The formulas for this are not terribly 
useful, but we can illustrate the results with a simple example. Let 
A^'(cp) = sin(aj) — sin(4> + aj), where j =1,2 for two oscillators. 
Then: 

0P i* = c \ T~i WT—2 ^T^f (25) 

1 + (t 2 + l) [sin 2 aj + sir a k \ 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



August 2013 | Volume 7 | Article 113 | 7 



Zhou et al. 



Correlated colored noise-induced synchrony 



Thus, for 0 < a\ < ci2 < n/2, we always have OPn > OP12 > 
OP22 for all x and sufficiently small values of c. This provides a 
simple and surprising illustration that heterogeneity will improve 
synchrony at low correlations for very simple PRCs. We remark, 
however, that this phenomenon does not always hold. Pairs of 
PRCs can be found such that OP is always bigger for both of 
the homogeneous oscillator pairs than for the heterogeneous 
oscillator pair, as can be seen from Figures 4F,G. 

3.2.3. PRC heterogeneity tunes the sharpness and peak position of 
the phase difference density 

If two neurons are identical but driven with partially correlated 
noise, then the peak of the phase difference density will be cen- 
tered at 4> = 0, which means that the two oscillators will tend 
to have the same phase. However, with heterogeneity, the peak 
will shift depending on the degree of heterogeneity, just as two 
coupled oscillators will shift if they have different intrinsic fre- 
quencies. Figure 5 shows how the peak of the phase difference 
density is shifted by oscillator heterogeneity. Using the two term 
double sinusoidal form PRC (Equation 63), we keep PRC 1 con- 
stant («i =0.1,£>i = 0.32) as we vary PRC 2 (&2 = 0.3 is constant 
and cii varies from — n to jt). From the results shown in Figure 5, 
we can conclude that heterogeneity can tune oscillator synchro- 
nization in both the sharpness and peak position of the phase 
difference density, which might be useful in neural signal coding. 
We also note that OP is minimized when the peak is at ±jt/2 and 
that "changing the sign" of the PRC (e.g., setting a 2 = it) shifts 
the peak but has very little effect on the OP. 

3.3. FREQUENCY DIFFERENCES HIGHLY LIMIT SYNCHRONIZATION 

In the above results, we assume that all oscillators have the 
same natural frequency, which means to = 0. This is a somewhat 
unreasonable assumption for real neurons. Thus, we now study 
how synchronization is dependent on the frequency differences 
between oscillators. Figure 6 shows the effects of frequency dif- 
ferences on a pair of oscillators that have different PRCs (of the 
two term double sinusoidal form, Equation 63) and are driven by 
partially correlated noise. With no frequency difference, the het- 
erogeneity in oscillator PRCs yields a shift in the peak position 
(Figure 5), consistent with previous measurements of synchrony 
between irregularly firing neurons (Tchumatchenko et al, 2010). 
This means that, if frequency differences can shift the peak in the 




-3-2-1 0 1 2 3 



FIGURE 5 I PRC heterogeneity tunes the sharpness and peak position 
of the phase difference density. OP vs. 82 (black axis) and the peak 
position of the phase difference density vs. (grey axis). Parameters used: 
a q = 0.1, b-, = 0.32. 6 2 = 0.3, x = 1, and c = 0.8. 



opposite direction [e.g., see Figure 1C of Burton et al. (2012)], 
then changes in frequency could "cancel" the effects of PRC het- 
erogeneity so that the peak of the phase difference density is at 0. 
This cancellation can be seen in Figure 6 near w = 0.2. However, 
this cancellation comes at a loss to precision, as seen by the 
decrease in OP. While not shown, we remark that the drop in 
OP is symmetric about oj = 0; thus, a negative frequency differ- 
ence will not result in a larger OP. While it remains to be proven, 
we conjecture that the OP is always maximal when there is no 
frequency difference. This differs from the case that we looked 
at in the previous section where heterogeneity (in PRCs) can 
sometimes lead to a larger OP than homogeneity. 

3.4. CORRELATED NOISE-INDUCED SYNCHRONIZATION IS 
DEPENDENT ON THE TIME CONSTANT OF THE NOISE 

Because of the natural decay times of synapses, broadband inputs 
into neurons have some temporal correlations. Thus, we now 
explore how the temporal properties of noise interact with het- 
erogeneities in the PRCs. Figure 7 shows that synchronization 
decreases monotonically as x increases for w = 0, while there 
exists an optimal value of x that achieves the greatest synchroniza- 
tion for w = 0.5. This means synchronization of two oscillators 
with different frequencies (i.e., co / 0) can have a resonance 
with x. Furthermore, as seen in Figures 7B,D the peak of the 
phase difference density depends on x only when there is a fre- 
quency difference between the two oscillators. In Figure 8, we 
explore this resonance in more detail where R($) is plotted as x 
varies. The left panels (showing the solution to the BVP and the 
results of Monte Carlo simulation) show that when w = 0, the 
peak position of R(<\>) is largely unchanged and the magnitude 
decreases monotonically with x. There is a sharp drop off in R(<\>) 
at x ~ 2. A different result emerges in the right panels, where a 
frequency difference exists (co = 0.5). At low and high values of 
x, R(<$>) is almost flat with a distinct resonance when x ~ 1. We 
see the same resonance in the biophysical ML model when the 
neurons have different frequencies and different PRCs (Figure 9). 

We can see why the frequency differences are needed for res- 
onance by considering the simplest example of identical PRCs 
of the form A((J>) = sin a — sin((f> + a). In this case, we solve the 



BVP: 




FIGURE 6 | Frequency differences limit noise-induced synchronization. 

OP decreases quickly as frequency differences increase (black axis). The 
peak position of the phase difference density is shifted by changing 
frequency differences between two oscillators (grey axis). Parameters used 
here: a, = 0.1, b, = 0.32, a 2 = 0.6, = 0.3, x = 1, and c = 0.8. 
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FIGURE 7 | Frequency differences between oscillators change the 
dependence of synchrony on the time constant of the correlated noise. 

(A) OP and (B) the peak position of the phase difference density for 
oscillators with no frequency difference (u> = 0). (C) OP and (D) the peak 
position of the phase difference density for oscillators with a frequency 
difference (to = 0.5). Parameters used: a-\ = 0.1, b\ = 0.32, a 2 = 0.6, 
= 0.3, and c = 0.8. 



FIGURE 9 | Frequency differences between biophysically realistic 
spiking neuron models change the dependence of synchrony on the 
time constant of the correlated noise. (A) OP and (B) peak position for 
model neurons with different PRCs but the same frequency 
(/, = 120, <]>! = 0.04. k = HO, and o> 2 = 0.04616). (C) OP and (D) peak 
position for model neurons with different PRCs (same as above) and slightly 
different frequencies (/, = 120, o>i = 0.041, l 2 = 110, and = 0.04616). 




FIGURE 8 | Synchronization of oscillators with different frequencies 
resonates with the time constant of the correlated noise. (A) BVP 

solution and (C) results of Monte Carlo simulation for f?(o» between two 
oscillators with no frequency difference (to = 0) as t varies. X-axis shows 
the phase difference, o), and higher fl(<])) is plotted as hotter colors in the 
heat map. (B) and (D), same as (A) and (C) for two oscillators with a 
frequency difference (to = 0.5). Parameters used: a-] = 0.1, b-j = 0.32, 
a 2 = 0.6, t>2 = 0.3, and c = 0.8. 



where G(4>, x) = (1 + x 2 ) sin(a) 2 + 1 — ccosc)). Here, a is pro- 
portional to the frequency difference. In particular, note that 
when oj = 0, G is independent of x and otherwise, x acts 
to weaken the correlated noise-induced synchronization as it 
increases the part of G that is not phase dependent. However, 



the right side of this equation shows that the effect of the fre- 
quency difference is minimized when x = 1, and thus we expect 
resonance in the OP. This effect disappears when a = 0. 

4. DISCUSSION 

In our current study, we have extended a number of previous 
results describing the ability of neural oscillators to synchronize 
in the presence of correlated noise. Our methods are similar to 
those in Burton et al. (2012), with the additional aspect that we 
now use colored noise (OU process). The properties of the noise 
show up only through a convolution of the autocorrelation func- 
tion of the noise with the phase functions h nm (4>) that, in turn, 
depend only on the PRCs (see Appendix, Equation 56). Thus, we 
could easily generalize this work to noise with an arbitrary auto- 
correlation function. In addition, we have now included many 
more examples of the theory and shown that the conclusions 
from the perturbation theory continue to be valid for full bio- 
physical models (cf. Figure 2). Further, we have shown that for 
low input correlations, heterogeneity can actually improve syn- 
chrony both pairwise and in large populations. We demonstrated 
that this theoretical effect can be seen in experimental recordings 
of regularly firing olfactory bulb mitral cells. Thus, we have sig- 
nificantly extended the findings presented in Burton et al. (2012), 
and our results on colored noise further suggest some experimen- 
tally testable phenomena, such as the resonance seen in slightly 
detuned oscillators (Figures 7-9). These novel findings and their 
biological implications are discussed in more detail below. 

4.1. HETEROGENEITY CAN IMPROVE SYNCHRONY 

We found that correlated noise can synchronize a heterogeneous 
pair of oscillators (comprised of a "bad synchronizer" and a 
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"good synchronizer") better than a homogeneous pair of bad 
synchronizers at low levels of input correlation and verified this 
experimentally. We showed that good (bad) synchronizers are 
characterized by having a relatively low (high) DC component in 
their PRC. Consistent with this, oscillators with the generic "type 
II" PRC (i.e., sine))) are better synchronizers than oscillators with 
the generic "type I" PRC (i.e., 1 — cosc|)). 

Several authors have previously studied the effects of het- 
erogeneity on the transfer of correlation. As we noted in 
Introduction, at low correlations, the output correlation is lin- 
early proportional to the input correlation through a factor, S, 
called the susceptibility (de la Rocha et al., 2007; Shea-Brown 
et al., 2008). If we let S(A, B) denote the susceptibility for two 
neurons, A, B, then what we have found in our current study 
is that in some cases, S(A,A) > S(A, B) > S(B,B). Note that 
in our study, we are looking at output correlation related to 
spike-to-spike synchronization, whereas in many other studies 
of output correlation, the interest is in spike count correlation. 
We can regard our measure of synchrony as the same as spike 
count correlation, but over a time window that is of the order 
of the mean interspike interval. In a recent paper, (Shea-Brown 
et al., 2008) showed that for spike count correlation, S(A, B) = 
JS(A,A)S(B,B) and thus, trivially, we can obtain S(A, A) > 
S(A, B) > S(B, B) when A is "better" than B at transferring cor- 
relation. We want to emphasize that their result is for long time 
windows (that is, the window length tends to infinity). Which 
neurons are better than others at the transfer of correlation 
depends very strongly on the window of time through which 
you measure the correlation. Indeed, Barreiro et al. (2010) and 
Abouzeid and Ermentrout (201 1) showed that type II PRCs have 
larger susceptibilities than type I for short time windows (i.e., 
spike synchrony) but the trend is reversed for large time windows 
(i.e., rate correlation). 

Interestingly, the efficiency of correlated-noise induced syn- 
chronization is also modulated by firing rate in the low input 
correlation regime (de la Rocha et al., 2007; Tchumatchenko 
et al., 2010). Given that changes in firing rate can modulate PRC 
shape in biophysically realistic neuron models and in real neurons 
(Gutkin et al, 2005; Marella and Ermentrout, 2008; Stiefel et al, 
2008, 2009; Schultheiss et al, 2010; Fink et al, 2011; Burton et al., 
2012), whether or not (and the degree to which) PRC hetero- 
geneity will enhance synchrony may depend in part on the firing 
rate. However, in the simplest cases (such as models like the leaky- 
integrate and fire model and the theta model), the only effect of 
the firing rate on the shape of the PRC is to change its amplitude. 
Since amplitude (but not shape) changes can be absorbed into the 
size of the noise, and our theory shows that the phase difference 
density is independent of the size of the noise (at least, if it is small 
enough), changes in firing rate will have no effect on the synchro- 
nization of pairs of neurons firing at the same or nearly the same 
rates. 

The ability of cellular heterogeneity to regulate which oscilla- 
tors synchronize best as a function of input correlation likely con- 
tributes to coding in many neural systems. In the olfactory bulb, 
where oscillatory synchrony appears to be critical to olfactory 
coding [for review, see Bathellier et al. (2010)], tens of "sister" 
mitral cells are linked to each glomerulus where they receive 



highly correlated afferent input (Carlson et al, 2000; Schoppa 
and Westbrook, 2001). Each sister mitral cell of a glomerulus may 
also participate in independent (i.e., unshared) lateral inhibitory 
circuits with non-sister mitral cells of surrounding glomeruli, 
mediated by local inhibitory granule cells (Dhawale et al., 2010; 
Tan et al., 2010). On average, sister mitral cells are thus subject to 
high input correlations while non-sister mitral cells are subject to 
low (though non-zero) input correlations (Dhawale et al., 2010). 
Further, we and others have demonstrated that mitral cells exhibit 
substantial cell-to-cell heterogeneity (Padmanabhan and Urban, 
2010; Angelo and Margrie, 201 1; Angelo et al, 2012; Burton et al, 
2012). Based on our current results, this heterogeneity will thus 
act to reduce output synchrony of sister mitral cells but enhance 
output synchrony of non-sister mitral cells. Thus, in the context 
of the olfactory system, heterogeneity will promote encoding of 
combinatorial sensory information (i.e., activation of non-sister 
mitral cells by odor combinations). 

Our results suggest that heterogeneity can only enhance 
correlation-induced synchronization by a moderate amount 
between two neural oscillators (up to 36% in BVP solutions 
using mitral cell PRCs). Two properties of the olfactory bulb 
nevertheless suggest that even this moderate effect can signifi- 
cantly influence patterns of oscillatory synchrony in the olfac- 
tory system. First, the reciprocal dendrodendritic connectivity 
between mitral cells and granule cells enables activity-dependent 
regulation of granule cell recruitment (Arevian et al., 2008), 
which can lead to amplification of granule cell-mediated corre- 
lated noise-induced synchronization (Marella and Ermentrout, 
2010). Second, mitral cells separated by up to ~2 mm can 
engage in lateral inhibitory interactions (Egger and Urban, 
2006), thus multiplying the synchrony-enhancing effect of cel- 
lular heterogeneity across a potentially large fraction of the 
~40,000 total mitral cells per mouse olfactory bulb (Benson 
et al., 1984). Whether neural oscillator heterogeneity exists in, 
and significantly enhances, correlated-noise induced synchrony 
in other brain regions remains a promising topic of future 
research. 

4.2. RESONANCE 

In addition to the above findings, we found that there exists 
some resonance of correlated noise-induced synchronization with 
respect to the time scale of the noise. That is, we found a local 
maximum in OP as the time scale of the correlated noise var- 
ied. Surprisingly, this only occurs when there is a difference in 
the frequencies between the two oscillators. The requirement for 
the frequency difference would seem to contradict earlier work 
(Galan et al., 2008), where it was found that the Liapunov expo- 
nent was most negative when the noise has a particular time scale. 
However, when the noise is only partially correlated, the uncor- 
rected part of the noise causes a drift in the phase difference. 
The degree of this drift is also dependent on the time scale of 
the noise, and thus the two effects cancel. A frequency difference 
breaks this symmetry by adding an additional drift term, which 
prevents one from factoring out the resonance. A frequency dif- 
ference thus leads to a dependence of OP on the time scale of 
the noise. We have not yet tested this idea experimentally, but it 
seems to be quite robust, having been found in both the simple 
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phase models (Figures 7, 8) and in the biophysical ML model 
(Figure 9). 

4.3. LIMITATIONS OF THE THEORY 

The analysis that we have done in this paper and in our earlier 
papers requires that the neurons fire almost periodically. This 
means that the activity of the neurons is mean driven rather than 
fluctuation driven so that the coefficient of variation of the inter- 
spike intervals should be small. While this may not be the case in 
all areas of the brain, there are many regions, such as the olfactory 
bulb, where the firing rate can be quite regular and synchronous 
as indicated by the large rhythmic local field potentials. Assuming 
that the neurons are firing at a fairly regular rate, it is also reason- 
able to ask how well the PRC describes such noisy neurons. An 
extensive review of the caveats of PRC theory for real neurons can 
be found in Smeal et al. (2010). Another issue is the actual esti- 
mation of the PRCs in the presence of noise. Several studies have 
shown that background synaptic activity and other forms of noise 
do not significantly affect the shape of the PRC (Ermentrout et al., 
2011; Netoffetal., 2012). 



In conclusion, we have extended our previous work 
to demonstrate that oscillator heterogeneity and frequency 
differences interact with the time scale of input noise 
to regulate how correlated noise synchronizes uncoupled 
oscillators. 

4.4. DATA SHARING 

All codes are available by request from the authors. They include 
Matlab and XPPaut files. 
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APPENDIX 

REDUCTION TO A PHASE MODEL 

Consider a general oscillator receiving a possibly noisy time-dependent signal: 

^ = F(X) + cN(X, f) (27) 
dt 

Here N(X, t) is the external or imposed inputs into the system. For single compartment neural models, N will typically only affect the 
membrane potential, e.g., as an injected or synaptic current. We assume that X' = F(X) has as a solution an exponentially stable limit 
cycle, U(t + X) = [/(f) and that € is a small positive parameter characterizing the magnitude of the input. We are interested in how 
the phase of the limit cycle evolves in time in the presence of small inputs. The phase of a limit cycle is easy to define when a point 
lies on the limit cycle. For example, for neurons, the phase is just the time since the last spike of the cell. However, if the limit cycle is 
attracting, then it is also possible to define the phase of a point that is near, but not directly on the limit cycle. Specifically, there is a 
function ©(X) that maps a point near the limit cycle, X, to the phase that it will eventually reach as it is attracted to the limit cycle 
(asymptotic phase). Clearly ©([/(f)) = t. Define the phase to be G(f) = 0(X(f)), so that by the chain rule: 

dQ dX 

-=V,8(X).- (28) 

= V X @(X) ■ F(X) + €V x 0(X) • N(X(t), t) (29) 

= l + tV x ®(X)-N(X(f),f) (30) 

Thus, in the absence of inputs, the phase moves around the circle at constant velocity. This expression is exact, but not very helpful since 
it requires knowledge of the solution X{t) . Kuramoto's approximation (which is valid for small €) is to replace X (t) in the right-hand 
side by [/(9(f)), where U is the unperturbed limit cycle (Kuramoto, 2003). This closes the system yielding: 

^- = l + eZ(Q)-N(U(Q),t) (31) 
dt 

where we have defined Z(9) := Vx® (U(Q)). The function, Z(9) is the so-called adjoint function satisfying the linear equation: 

Z' = - (D x F(U(t))f Z(t) (32) 

with Z(t) ■ U'(t) = 1. Here DxF(U(t)) means the linearization ofF(X) evaluated along the limit cycle. 

In single compartment neuron models, inputs appear only in the voltage component of the neural oscillator in the form of external 
currents so that the dot product in Equation 31 becomes scalar multiplication: 

^ = 1 +eA(Q)I(U(Q),t)/C (33) 
dt 

where I is the input current, C is the capacitance, and A (6) is the voltage component of the vector Z. The quantity, A(G), is sometimes 
called the infinitesimal PRC and, for small perturbations, is proportional to the PRC. 

DERIVATION OF THE STATIONARY DENSITY OF PHASE DIFFERENCES 

The Langevin equations that drive the phase models (Equations 4-6) correspond to a forward Fokker-Planck (FP) equation that can 
be written as (Risken, 1984): 



3p 113/13 \ 3/13 \ 3 2 
1 1 +x) + — [ -— +y) +c- 



3f x [ dx \ 2 dx J dy \2 dy J dxdy 



p P (34) 

H 39 



id d 

-e — [Ai(9)*p] + rr[(A 2 (e + ®y - Ai(9)x)p] 
I 39 3q> 

2 9 
—€03 p 

When the distribution of phase differences is stationary, — = 0. Our goal is to exploit the smallness of € to compute this stationary 

3f 

density with which we can compute the marginal distribution of the phase difference. 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



August 2013 | Volume 7 | Article 113 | 13 



Zhou et al. 



Correlated colored noise-induced synchrony 



ANALYTICAL SOLUTION 

We expand the steady state p in €: 



We define the operator: 



,3P 



p(x,y, 0, (t)) = po(x,y, 0, to + api(x,y, 0, to + t 2 p 2 {x,y, 0, to 

fill P0 ^ x,y ' 9 ' ^ dxd y dQd § - 1 
//// pn(x ' y ' e ' ^^y^^ = °' 



(35) 



« = 1,2 



ira/i 



3x 



+ x + 



3/13 



3y \2 3y 



3x3y 



+ 



30 



(36) 



At steady state condition ( — =0), we substitute the above expansion into the FP equation and collect the powers of €. We need to go 
3t 



0 = LnPo 
0 = L 0 Pi 

0 = L 0 p 2 



9 3 

— [Ai (0)*po] + — [( A 2 (9 + +)y - A , (9)*) p 0 ] 
30 3® 



30 



[Ai(9)xpi] + — [(A 2 (9 + toy ~ Ai(0)*)pi] - a—p 0 



Solving Equation 37 

Equation 37 is just a linear separable equation, independent of 4>, so, by inspection: 



(37) 
(38) 

(39) 



p 0 (x,y,Q,to = —G(x,y)R(to 
lit 



where: 



G(x, y) 



-^(x 2 +y 2 -2cxy) 



(40) 



(41) 



and R(to remains to be determined. Note that G(x, y) is just the standard stationary solution to the multivariate OU equation. At 
this juncture, we remark that our main goal is to find -R((|)), which is the marginal density of the phase differences between the two 
oscillators. 

Solving Equation 38 

Both Equations 38 and 39 have the form LoP = b(x, y, 0). in operates on the space of functions defined on R 2 x S 1 that are twice 
continuously differentiable in x, y and continuously differentiable in 0. In this space, Lq has a one-dimensional nullspace spanned by 
G{x,y) (constant in 0) and so Lq is not invertible. However, Lop(x, y, 9) = b(x, y, 9) does have a solution provided that b(x, y, 9) is 
orthogonal to the null space of Lq, which is the adjoint linear operator of in. Since In is a standard probability operator, its adjoint is 
always 1 (i.e., the function that is 1 everywhere). 



Since: 



h(x, y, 0) = * [ X ' y) [A[(B)R(.to ~ Ai(e)R'((j))] + yG(x, y)[A' 2 (Q + toR(to + A 2 (9 + 
lit 



we see that JJJ b\ (x, y, Q)dxdydQ = 0. Thus, In Pi = b\ has a solution. Since: 

Lq[xG(x, y)] = —xG(x,y)/x 
L 0 [yG(x,y)] = -yG{x,y)/x 



we look for a solution of the form: 



Pi(x, y, 0, to 



W\ (9, cj))xG(x, y) + W2(9, toyG(x, y) 
2ji 



(42) 

(43) 
(44) 

(45) 
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Inserting this into Equation 38, we find that wi(Q, fy) must satisfy: 

3 Wi(0,cb) 
at) x 

4w 2 (6, (» + W2(6 '^ = - A' 2 (9 + (t.)i?((t)) - A 2 (9 + c|>)i?'($) (47) 

Wj must be periodic functions of 6; we defer their exact solution to later, but note that there is always a unique periodic solution to 
each of these equations. 



Solving Equation 39 

We now have: 

9 9 9 

b 2 (x, y, 9) = — [A^xpi] + — [(A 2 (9 + <\>)y - Ai(9)x)pi] + a— p 0 (48) 
99 9(p 9(p 

In order to solve Equation 39, for p 2 (x, y, 9, (()), we must have: 
0= / h(x,y, Q)dxdydQ 

J J-oo JO 

9 f /•/• 00 f 2 ' 1 f A 2 (9 + (b) , 

Ai(9) , a 1 

—[wi(Q, (t))x 2 + w 2 (9, <$>)xy]G(x, y) H _R((t))G(x,y) [ dxdydQ 

2tt 2tt J 

1 9 



2it 

{A 2 (9 + ct>)[w 2 (9, (()) + c • wi(9, ct>)] - Ai(9)[wi(9, 4>) + c ■ w 2 (9, (t>)]}d9 + 4nu>R(<b) 



o 



4tt 9(() 

^-^-[f(<\>) + 4naRm (49) 
4jt 9A 



where: 



/(<t>)=/ [A 2 (9+^(9,(|))-A 1 (9)v 1 (9,<t))]d9 (50) 
Jo 

Vl (9, (t>) = Wi(9, (t>) + cw 2 (9, (t>) 
v 2 (9, (t>) = w 2 (9, (t>) + cwi(9, (|)) 

Given Equations 46 and 47, we see that v\ (9) and v 2 (9) satisfy: 

Vi(9)+ — = -[A;(9)4- c . A' 2 (9 + (())]R(c|)) + [Ai(9)-c- A 2 (9 + (|))]R'((|)) 

:=«i(6) (51) 

V 2 (Q) + ^ = -[c ■ A; (9) + A^(9 + cJ>)]R(d>) + [c ■ Ai(9) - A 2 (9 + c|>)]R'(<|>) 

:= ?2(9) (52) 

For Equations 51 and 52, we can write down the solution of V\ (9) and v 2 (9) in terms of q\ (9) and q 2 (9) (see Appendix): 



v„(9) = / e~^q n (6- s)ds, n=l,2 (53) 
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We substitute v„(4>) into f(<\>), 



/(<))) = / [A 2 (G + (|))v 2 (G) - Ai(9)vi(6)]d9 
Jo 



poo p2it 

/ e -Us [A 2 (e + <|))52(e-s)-Ai(9)^(e-s)]de 

Jo Jo 



where: 



poo pin 

= e~Us h(Q,<b,s)dQ (54) 
Jo Jo 



h(Q, ct>, s) = A 2 (9 + $)q 2 (B - s) - Ai(9)g 1 (9 - s) 
= Ai (6 - s)[c • A 2 (G + *) - Ai(9)]R'(c|)) 

+ A 2 (9 + * - s)[-A 2 (9 + 4>) + c • Ai(9)]i?'(<|>) 
+ A; (9 - s)[-c • A 2 (9 + <j)) + Ai(9)]_R(())) 

+ A' 2 (9 + c|) - s)[-A 2 (9 + (()) + c • Ai(9)]*(4>) (55) 



Define: 



vmn(«l>) = / h mn (s + §)e~ids (56) 



1 



fe m „(s) = / A m (9)A„(9 + s)d9 
Jo 



Since Ai(9) and A 2 (9) are periodic functions, 



where: 



rZn pin 

/ ft(9. 4>, s)dB = / {A 1 (9-5){[cA 2 (9 + c|))-Ai(9)]R'(c|)) + [cA^(9 + (|))-A' 1 (9)]J?(<|))} 
Jo Jo 

+ A 2 (9 + (|> - s){[cAi(6) - A 2 (9 + ct>)]£'(cf>) - [cA'^G) - A' 2 (9 + cp)] J R(cp) }}ciG 
= {c[hn(s + <|>) + h 2 i(s - +)] - [/Jii(s) + felCs)]}*^) 

+ {^feC 5 + +) + h 2l( S ~ +)] " ^ n (s + +) - ^ + *)]| + = 0 J *(4>) (57) 

/(((,) = / e- 'xds / fr(9, 4>, s)rf9 
Jo Jo 

= {cfeia C<t>) + «zi (— 4»)] " feii(0) + g 2 2(0)]}i?'(<|>) 

+ jc- ^feu(40+fci(-4>)] - l/n(0)-/ 22 (0)]}R(« 

= ^{[c • *(«>) - Ci]K(«|))} - C 2 R(® (58) 



(-♦) 09) 
Ci = gn(0)+g 22 (0) (60) 
C 2 = g / u (0)-g 22 (0) (61) 
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Combined with Equations 49-58, we have a boundary value problem (BVP): 

d 



,Alc ■ *(<M " Ci]K(4>)} + (47tco - C 2 )R(4>) = K (62) 



r 



J?((|)) = R(<\> + 2jt) 
*(4>) = + 2«) 
l?((t))# = 1 

C 2 

K = 2co - 

2 it 



To solve this BVP, we need to compute g((|>) for a given PRC. We use two forms of the PRC in this paper: 

A (9) = sin (a) - sin(6 + a) + b sin(29) (63) 

and 

A(9) = A[sin(5) - sin(9 + B)]e c(9 ~ 2lt) (G e (0, 2n), A (9) = A(9 + 2it)) (64) 

The required integrals can be computed for both PRC forms. More generally, all PRCs can be written in Fourier form and, again, the 
integrals are readily computed to obtain g(<\>) (see below). 

Small correlation approximation for /?(<(>) 

We use a BVP solver to get the numerical solution for R(<\>), but we would like to better understand the form of R(<\>) at low values 
of correlation, c, so we expand _R as a series in c. As K is dependent on c, we must also expand K in c. Finally, we need to keep the 
normalization condition for R(§), hence: 

R($) =R 0 (<)>) + d?i((|)) + ... K = K 0 + cK 1 + ... (65) 

R 0 ((|))=R(<|))| C=0 , f i?o(*)#=l 
J— it 



R n (§) = 0, n > 1 



We substitute these expressions into the BVP, Equation 62 and find after equating powers of c: 

- di^Ccb) + (4jtco - Ca)Ko(<» = *5> (66) 
-Ci^(^) + (4™ - C 2 )Ri(<|)) + feWRoC*)]' = *i (67) 

We can integrate both left and right side over [—it, jt], then use the assumptions and periodicity requirements above to get: 



Rewriting these equations, 



4ttw — C? , , 

= -, Ki = 0 (68) 

2n 



R' 0 (<\,)+DR 0 (<b) = Qo (69) 

Co — 4jtW Kn 

D=— , Q 0 = 

R'i($) + DRi(fy) = Qi (<(>) (70) 



Qi(*) 



C, 



we see: 



2jt 

[Ri(^)e D *]' = Qi((We z 
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We can use numerical methods to get the solution of i?i (c{)) given -Ro(^) and for some choices of PRCs we can get exact expressions. 
For example, for the two term double sinusoidal form PRC (Equation 63) we get: 



Ki(4>) 



"I — 

Ci 1 x 2 + 1 



cos(<\> + a.2 — «i) — Dsin(4> + «2 — «i) 2b\b2X 2 cos(2(|>) — Dsin(2c(>) 



l + D 2 



+ 



Ax 2 + 1 



4 + D 2 



(71) 



DIFFERENT PRCS 
Double sines form PRC 

For the PRC: 



We have: 



A m (9) = sin(a m ) — sin(G + a m ) + b m sin(26) 



h mn (s) = f A m (0)A„(9 + s)d9 
Jo 

= 2 it sin(a m ) sin(a„) + 7t cos(s — a m + a n ) + b m b n ix. cos(2s) 



POO 

,(<|>) = / h mn (s + $)e~*ds 
Jo 

= lux sin(a m ) sin(a„) + ttt. 



cos((() + a n — a m ) — x sin(cj) + a n — a m ) 



x 2 + 1 



+ b m b n Tix 



cos(2c|)) — 2x s'm(2<\>) 
4x 2 + 1 



, x cos(<\> + a n — a m ) + sin((() + a„ — a m ) 2x cos(2())) + sin(2())) 

?mnW = ~ 7Tt TT 2b m b„itx- 



x z + 1 



4x 2 + 1 



g(<W = guW +g2l(-$) 



cos(c|) + U2 — a\) cos(2(|)) 
4jtt, sin(«i) sin(«2) + 2nx ; h IbibjTiT- 



x 2 + l 



' 4x 2 + 1 



?\<\>) = g[ 2 (4>) ~ gn(-& 

sin((() + «2 — «i) 



-2ltT- 



T 2 + 1 



4b\b2Ttx 



sin(2(|)) 
4t 2 + 1 



Ci =^ii(0)+ & 2(0) = 2Ttt(sin 2 (fl 1 ) + sin 2 (fl 2 ))+^— + 1 



x 2 + 1 4t 2 + 1 



47tt 2 (fo 2 - fo 2 ) 

C2=/n(0)-^ 2 (0)= 4 y +1 V 



(72) 



(73) 
(74) 

(75) 

(76) 
(77) 
(78) 



Exponential-sine form PRC 

For empirical PRCs: 



We have: 



Ai(G) = fli[sin(fo 1 ) -sin(foi +6)]e Cl(9 - 27t) 
A 2 (9) = fl 2 [sin(fo 2 ) - sin(fo 2 + 0)]e C2(9 - 27t) 
9 e (0, 2:t) 



h mn (s) = Bi ■ e- c " s (hc u fcvi(s)) + B 0 • e c " s (fcc 0 , hv 0 (s)) 
gm»(4>) = Co • - e- c ""*(gc 1; gvi («>)) - e c »% 0 , gv 0 (*)) 



(79) 
(80) 
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Where s e [0, 2n) and (|> G [0, 2tc), both have the period of 2n. Note also that (,) means the inner product of two vectors. The above 
quantities are defined as: 

B\ = a m ■ a n ■ e iJtc " • [1 

B 0 = a m ■ a„ • [1 - e~ 2llCn ] 
1 

Di = c m + - 
x 

1 

Do = c n 

x 

sin(Z? m ) • sin(fc„) sin(Z? m ) . 

fc 0 = — — - — - — r^LCcm + c„) • sm(b n ) - cos(fo„)] 

, 1 , 1 , c m -\- C n 

h = 



2(c m + c„)' 4 + (c m + c n ) 2 ' 2[4 + (c m + c„) 2 ] 

(c m + c n ) ■ sin(fo„) sin(fo„) 
k 4 = , . ; ; ^— . h 



1 + (c m + c„) 2 ' ' 1 + (c m + c„) 2 
sin(fo m ) • sin(fo„) sin(b„) 

c m -\~ c n 1 -f- (c m + C M ) 2 



1 1 Cfi 



n = w} — ; — \ ' n = --r—, — ; — z?, n 



2(c m + c n )' 4 + (c m + c n ) 2 ' 2[4 + (c m + c„) 2 ] 

(c m + c n ) ■ sm(b m ) . sm(b m ) 



Ji 1 + (Cm + c„) 2 ' 15 1 + (c m + C„) 2 

hc x = [ko, h, k 2 , k } , k A , fc 5 ] 

hc 0 = [jo,ji,h,hJi,j5] 

hv\ = [1, cos(s + b n — b m ), sin(s — b m — b n ), cos(s — b m — b n ), sin(s — b m ), cos(s — b m )] 
hvo = [1, cos(s + b„ — b m ), sin(5 + fo,„ + cos(s + b m + b n ), sin(s + £>„), cos(s + b„)] 
Bi 1 + D? 

gci = -=■[ — -k 0 , -D x ki, ky, k3 - Difa, -fc 2 - Di*3, fc5 - Difc*, — fc4 - Di* 5 ] 

l+D\ Di 

Bo 1 + D 2 

gQ> = -r [— =— V Oq7i,;i,;3 + D 0 ; 2 , -72 + A>;3,;5 + D 0 ; 4 , -;4 + D 0 ; 5 ] 

I • />;; D 0 

gvi($) = [l, cos(4) + b n - b m ), sin((|) + fo„ - fo m ), sin((|) - b m - b n ), 

cos(4> — b m — b n ), sin(4> — b m ), cos(s — b m )] 
gvoify) = [h cos(<|> + fc„ - b m ), sin((() + fc„ - fc m ), sin(<\> + b m + b„), 

cos((|) + b m + b n ), sin((|> + &„), cos(s + b„)] 

Co = ^[(e- 2ltCm - 1) • (gcugvm) + (e 2 * c " - 1) • (^ 0 ^v 0 (0))] 

1 - e^T" 



Fourier form PRC 

For the Fourier form of the PRC: 



A m (9) = a„ 

k= — oo 



/z m „(s) = / A m (9)A„(6)ds 
Jo 

2u 



rill 

ki , fa 
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7. a„ 



,—iks 



2tc } a„ y -*e 

k= — oo 
roc 

%mnW) = / h mn (s+ $)e~*ds 
Jo 

oo c 
= 2lt >, fl m, Jc«n, -jfc / 

, — Jo 

k= — oo 



(81) 



2jt ^ fl ^ fcfl "-- t (i : -ik)e- ik * 



k= — oo 



/c2 + JL t 



(82) 
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